CD4-Benchmark
Thomas W. Battaglia
2022-09-25
Introduction
Reads were first de-multiplexed with ea-utils (fastq-multx) with an allowed barcode mismatch of 1 (Default). De-multiplexed reads were then subjected to an alignment against the UniVec vector contamination database using the adapter removal tool within ea-utils (fastq-mcf). Barcodes were extracted by remove the last 12 nucleotides from each read.
Next reads were subjected to an alignment against the Oligo library using a DNA aligner. We used Bowtie2…
Data import & preprocessing
We are going to first import the merged counts tables that is outputted from the Nextflow pipeline. This table already has the raw counts summarised from the alignments and merged into a single matrix.
Import data
# Import CD4 counts
counts = readRDS('../data/CD4-counts.rda')
head(counts)## cd4_bcell_alone_1 cd4_cd8_mansc1_1 cd4_cd8_mansc1_2 cd4_cd8_mansc1_3
## oligo_1 9153 45 53 38
## oligo_2 15918 76 94 80
## oligo_3 7795 28 29 34
## oligo_4 9068 132 39 44
## oligo_5 8207 25 55 41
## oligo_6 5805 23 42 32
## cd4_cd8_mock_1 cd4_cd8_mock_2 cd4_cd8_snord73_1 cd4_cd8_snord73_2
## oligo_1 39 55 25 47
## oligo_2 88 88 35 108
## oligo_3 25 26 14 42
## oligo_4 46 36 23 43
## oligo_5 45 42 19 54
## oligo_6 35 35 11 37
## cd4_library_maxi_1 cd4_library_maxi_2 cd4_mansc1_1 cd4_mansc1_2
## oligo_1 1969 1955 5790 5149
## oligo_2 3357 2816 10287 8929
## oligo_3 1506 1595 4868 4243
## oligo_4 1604 1761 4565 3622
## oligo_5 1226 1334 4701 3600
## oligo_6 1251 1191 3228 3239
## cd4_mansc1_3 cd4_mock_1 cd4_mock_2 cd4_snord73_1 cd4_snord73_2
## oligo_1 5117 6086 1844 6762 7154
## oligo_2 9485 10287 2803 13097 12108
## oligo_3 4508 5181 1422 6088 5464
## oligo_4 4399 5286 1396 5650 6561
## oligo_5 3825 4086 1275 5025 4895
## oligo_6 3035 3849 1305 4501 3700
## cd4_snord73_3
## oligo_1 7346
## oligo_2 13418
## oligo_3 5936
## oligo_4 6067
## oligo_5 5383
## oligo_6 4226# Import CD4 metadata
metadata = readRDS('../data/CD4-metadata.rda')
head(metadata)## Sample Group Barcode
## cd4_library_maxi_1 cd4_library_maxi_1 DNA GGAATTC
## cd4_library_maxi_2 cd4_library_maxi_2 DNA TTAGTGG
## cd4_bcell_alone_1 cd4_bcell_alone_1 Bcell ACAACGA
## cd4_mock_1 cd4_mock_1 Mock ACCAATG
## cd4_mock_2 cd4_mock_2 Mock CAACAGC
## cd4_snord73_1 cd4_snord73_1 SNORD CACACGT# Import MultiQC results
report = readRDS('../data/CD4-multiqc.rda')
head(report)## # A tibble: 6 × 47
## Sample raw_total_seque… filtered_sequen… sequences is_sorted `1st_fragments`
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 cd4_bce… 32649199 0 32649199 1 32649199
## 2 cd4_cd8… 204206 0 204206 1 204206
## 3 cd4_cd8… 196064 0 196064 1 196064
## 4 cd4_cd8… 200980 0 200980 1 200980
## 5 cd4_cd8… 209112 0 209112 1 209112
## 6 cd4_cd8… 226005 0 226005 1 226005
## # … with 41 more variables: last_fragments <dbl>, reads_mapped <dbl>,
## # reads_mapped_and_paired <dbl>, reads_unmapped <dbl>,
## # reads_properly_paired <dbl>, reads_paired <dbl>, reads_duplicated <dbl>,
## # reads_MQ0 <dbl>, reads_QC_failed <dbl>, `non-primary_alignments` <dbl>,
## # total_length <dbl>, total_first_fragment_length <dbl>,
## # total_last_fragment_length <dbl>, bases_mapped <dbl>,
## # `bases_mapped_(cigar)` <dbl>, bases_trimmed <dbl>, …# Import oligo annotations
annotations = readRDS('../data/CD4-annotations.rda')
head(annotations)## # A tibble: 6 × 3
## oligo_id gene_name gene_id
## <chr> <chr> <chr>
## 1 oligo_1 MAGEA1 MAGEA1 #1
## 2 oligo_2 <NA> NA #1
## 3 oligo_3 <NA> NA #2
## 4 oligo_4 <NA> NA #3
## 5 oligo_5 <NA> NA #4
## 6 oligo_6 <NA> NA #5Mapping quality
# Mean
mean(report$raw_total_sequences)## [1] 11596529median(report$raw_total_sequences)## [1] 7196006range(report$raw_total_sequences)## [1] 114571 32649199Number of mapped reads
p1 = report %>%
select(Sample, reads_unmapped, reads_mapped) %>%
mutate(Sample = str_replace(Sample, '-bamstats', '')) %>%
pivot_longer(cols = c("reads_unmapped", "reads_mapped"), names_to = "feature", values_to = "reads") %>%
mutate(feature = if_else(feature == "reads_unmapped", "Unmapped", feature)) %>%
mutate(feature = if_else(feature == "reads_mapped", "Mapped", feature)) %>%
ggplot(aes(x = reorder(Sample, reads), y = log10(reads), fill = feature)) +
geom_col() +
coord_flip() +
theme_classic(base_size = 10) +
theme(legend.position = "top") +
scale_fill_brewer(palette = "Set2") +
xlab("") +
ylab("Number of reads (log10)") +
ggtitle("CD4: mapping rate")
p1
### Percent of mapped reads
p2 = report %>%
select(Sample, reads_unmapped_percent, reads_mapped_percent ) %>%
mutate(Sample = str_replace(Sample, '-bamstats', '')) %>%
pivot_longer(cols = c("reads_unmapped_percent", "reads_mapped_percent"), names_to = "feature", values_to = "reads") %>%
mutate(reads = reads / 100) %>%
mutate(feature = if_else(feature == "reads_unmapped_percent", "Unmapped", feature)) %>%
mutate(feature = if_else(feature == "reads_mapped_percent", "Mapped", feature)) %>%
ggplot(aes(x = reorder(Sample, reads), y = reads, fill = feature)) +
geom_col() +
coord_flip() +
theme_classic(base_size = 10) +
theme(legend.position = "top") +
scale_fill_brewer(palette = "Set2") +
scale_y_continuous(labels = scales::percent) +
xlab("") +
ylab("Percent of reads")
p2
# Plot together
p = (p1 | p2) + plot_layout(guides = "collect") & theme(legend.position = "bottom")
p
ggsave("figures/Benchmark-CD4-mapping-qc.pdf", p, width = 8, height = 5)Normalize data
# Match the sample ID's between samples
counts = counts[,match(metadata$Sample, colnames(counts))]
# Make DESeq2 object
ddsMat <- DESeqDataSetFromMatrix(countData = counts,
colData = metadata,
design = ~Group)
# Get normalize counts
ddsMat <- estimateSizeFactors(ddsMat)
# Filter flow counts
ddsMat.fil = ddsMat %>%
filter_low(percentile = 0.05)
# Get normalized counts
counts.norm.fil <- (counts(ddsMat.fil, normalized = TRUE) + 0.5) %>%
as.data.frame() %>%
rownames_to_column("oligo_id")
# Create a long table
counts.norm.fil.long = counts.norm.fil %>%
pivot_longer(cols = -oligo_id, names_to = "sampleId", values_to = "norm")Mock replicates
counts.norm.fil %>%
ggplot(aes(x = log10(cd4_mock_1), y = log10(cd4_mock_2))) +
geom_point(colour = alpha("grey", 0.7)) +
xlab("Mock replicate #1 (log10)") +
ylab("Mock replicate #2 (log10)") +
theme_light() +
stat_cor(method = "pearson", aes(label = ..r.label..)) +
geom_smooth(method = "lm", se = FALSE)## `geom_smooth()` using formula 'y ~ x'
### SNORD73
replicates
counts.norm.fil %>%
ggplot(aes(x = log10(cd4_snord73_1), y = log10(cd4_snord73_2))) +
geom_point(colour = alpha("grey", 0.7)) +
xlab("SNORD73 replicate #1 (log10)") +
ylab("SNORD73 replicate #2 (log10)") +
theme_light(base_size = 11) +
stat_cor(method = "pearson", aes(label = ..r.label..)) +
geom_smooth(method = "lm", se = FALSE)## `geom_smooth()` using formula 'y ~ x'
counts.norm.fil %>%
ggplot(aes(x = log10(cd4_snord73_2), y = log10(cd4_snord73_3))) +
geom_point(colour = alpha("grey", 0.7)) +
xlab("SNORD73 replicate #1 (log10)") +
ylab("SNORD73 replicate #2 (log10)") +
theme_light(base_size = 11) +
stat_cor(method = "pearson", aes(label = ..r.label..)) +
geom_smooth(method = "lm", se = FALSE) ## `geom_smooth()` using formula 'y ~ x'
counts.norm.fil %>%
ggplot(aes(x = log10(cd4_snord73_1), y = log10(cd4_snord73_2))) +
geom_point(colour = alpha("grey", 0.7)) +
xlab("SNORD73 replicate #1 (log10)") +
ylab("SNORD73 replicate #2 (log10)") +
theme_light(base_size = 11) +
stat_cor(method = "pearson", aes(label = ..r.label..)) +
geom_smooth(method = "lm", se = FALSE)## `geom_smooth()` using formula 'y ~ x'
MANSC1 replicates
counts.norm.fil %>%
ggplot(aes(x = log10(cd4_mansc1_1), y = log10(cd4_mansc1_2))) +
geom_point(colour = alpha("grey", 0.7)) +
xlab("MANSC1 replicate #1 (log10)") +
ylab("MANSC1 replicate #2 (log10)") +
theme_light(base_size = 11) +
stat_cor(method = "pearson", aes(label = ..r.label..)) +
geom_smooth(method = "lm", se = FALSE) ## `geom_smooth()` using formula 'y ~ x'
counts.norm.fil %>%
ggplot(aes(x = log10(cd4_mansc1_2), y = log10(cd4_mansc1_3))) +
geom_point(colour = alpha("grey", 0.7)) +
xlab("MANSC1 replicate #2 (log10)") +
ylab("MANSC1 replicate #3 (log10)") +
theme_light(base_size = 11) +
stat_cor(method = "pearson", aes(label = ..r.label..)) +
geom_smooth(method = "lm", se = FALSE) ## `geom_smooth()` using formula 'y ~ x'
counts.norm.fil %>%
ggplot(aes(x = log10(cd4_mansc1_1), y = log10(cd4_mansc1_3))) +
geom_point(colour = alpha("grey", 0.7)) +
xlab("MANSC1 replicate #1 (log10)") +
ylab("MANSC1 replicate #3 (log10)") +
theme_light(base_size = 11) +
stat_cor(method = "pearson", aes(label = ..r.label..)) +
geom_smooth(method = "lm", se = FALSE) ## `geom_smooth()` using formula 'y ~ x'
Plot normalized histogram
counts.norm.fil.long %>%
ggplot(aes(x = log10(norm))) +
geom_histogram(aes(y=..density..), binwidth=0.25, colour="black", fill="white") +
geom_density(alpha = 0.2, fill = "#FF6666") +
facet_wrap(sampleId~.) +
theme_light(base_size = 11) +
ylab("Density") +
xlab("Normalized reads per oligo (log10)") 
Determine dropout oligo’s
SNORD73 vs. Mock
snord_mock.res = compute_stats(ddsMat.fil, trt = "SNORD", ref = "Mock", fitType = "parametric", lfcThreshold = 0.25, alpha = 0.05) %>%
left_join(annotations) %>%
mutate(log2FoldChange = round(log2FoldChange, 2)) %>%
mutate(log2FoldChangeShrunken = round(log2FoldChangeShrunken, 2)) %>%
mutate(gene_name = if_else(gene_name == "RPS3Amut", "SNORDmut", gene_name)) %>%
mutate(gene_name = if_else(gene_name == "RPS3Awt", "SNORDwt", gene_name))## using pre-existing size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## -- note: fitType='parametric', but the dispersion trend was not well captured by the
## function: y = a/x + b, and a local regression fit was automatically substituted.
## specify fitType='local' or 'mean' to avoid this message next time.## final dispersion estimates## fitting model and testing## using 'normal' for LFC shrinkage, the Normal prior from Love et al (2014).
##
## Note that type='apeglm' and type='ashr' have shown to have less bias than type='normal'.
## See ?lfcShrink for more details on shrinkage type, and the DESeq2 vignette.
## Reference: https://doi.org/10.1093/bioinformatics/bty895##
## out of 4525 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.25 (up) : 0, 0%
## LFC < -0.25 (down) : 2, 0.044%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 14)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results## Joining, by = "oligo_id"# Significant oligos
snord_mock.res %>%
filter(padj < 0.05)## oligo_id baseMean log2FoldChange lfcSE stat pvalue
## 1 oligo_1889 195.6394 -2.50 0.2088072 -10.762587 2.584900e-27
## 2 oligo_4271 162.4496 -2.35 0.2284288 -9.182969 2.096831e-20
## padj log2FoldChangeShrunken gene_name gene_id
## 1 1.169667e-23 -0.83 SNORDmut RPS3Amut #1
## 2 4.744080e-17 -0.69 SNORDmut RPS3Amut #2# Browse
snord_mock.res %>%
select(gene_id, log2FoldChange, padj, gene_name, gene_id) %>%
DT::datatable(extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
))# Plot MA-plot
snord_mock.maplot = snord_mock.res %>%
ggplot(aes(x = log2(baseMean), y = log2FoldChange, label = gene_name)) +
geom_point(colour = alpha("grey", 0.7)) +
geom_point(colour = alpha("grey", 0.7)) +
geom_point(data = filter(snord_mock.res, gene_name %in% c("SNORDmut")), size = 3, colour = "#efc383") +
geom_point(data = filter(snord_mock.res, gene_name %in% c("SNORDwt")), size = 3, colour = "#a6a6a6") +
geom_text_repel(data = filter(snord_mock.res, gene_name %in% c("SNORDmut", "SNORDwt")), box.padding = 0.5, max.overlaps = Inf, size = 2.5) +
theme_classic() +
theme_classic() +
xlab("Mean counts (log2)") +
ylab("Fold change log2(SNORD73/Mock)") +
ggtitle("SNORD")
ggplotly(snord_mock.maplot)MANSC1 vs. Mock
mansc_mock.res = compute_stats(ddsMat.fil, trt = "MANSC", ref = "Mock", fitType = "parametric", lfcThreshold = 0.25, alpha = 0.05) %>%
left_join(annotations) %>%
mutate(log2FoldChange = round(log2FoldChange, 2)) %>%
mutate(log2FoldChangeShrunken = round(log2FoldChangeShrunken, 2)) %>%
mutate(gene_name = if_else(gene_name == "RPS3Amut", "SNORDmut", gene_name)) %>%
mutate(gene_name = if_else(gene_name == "RPS3Awt", "SNORDwt", gene_name))## using pre-existing size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## -- note: fitType='parametric', but the dispersion trend was not well captured by the
## function: y = a/x + b, and a local regression fit was automatically substituted.
## specify fitType='local' or 'mean' to avoid this message next time.## final dispersion estimates## fitting model and testing## using 'normal' for LFC shrinkage, the Normal prior from Love et al (2014).
##
## Note that type='apeglm' and type='ashr' have shown to have less bias than type='normal'.
## See ?lfcShrink for more details on shrinkage type, and the DESeq2 vignette.
## Reference: https://doi.org/10.1093/bioinformatics/bty895##
## out of 4525 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.25 (up) : 0, 0%
## LFC < -0.25 (down) : 2, 0.044%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 10)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results## Joining, by = "oligo_id"# Significant oligos
mansc_mock.res %>%
filter(padj < 0.05) ## oligo_id baseMean log2FoldChange lfcSE stat pvalue
## 1 oligo_1895 618.9236 -1.94 0.1207156 -13.993815 8.502020e-45
## 2 oligo_4277 429.8109 -2.05 0.3513332 -5.119629 1.530689e-07
## padj log2FoldChangeShrunken gene_name gene_id
## 1 3.847164e-41 -1.21 MANSC1mut MANSC1mut #1
## 2 3.463184e-04 -0.33 MANSC1mut MANSC1mut #2# Browse
mansc_mock.res %>%
select(gene_id, log2FoldChange, padj, gene_name, gene_id) %>%
DT::datatable(extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
))# Plot MA-plot
mansc_mock.maplot = mansc_mock.res %>%
ggplot(aes(x = log2(baseMean), y = log2FoldChange, label = gene_name)) +
geom_point(colour = alpha("grey", 0.7)) +
geom_point(colour = alpha("grey", 0.7)) +
geom_point(data = filter(mansc_mock.res, gene_name %in% c("MANSC1mut")), size = 3, colour = "#aec3a0") +
geom_point(data = filter(mansc_mock.res, gene_name %in% c("MANSC1wt")), size = 3, colour = "#a6a6a6") +
geom_text_repel(data = filter(mansc_mock.res, gene_name %in% c("MANSC1mut", "MANSC1wt")), box.padding = 0.5, max.overlaps = Inf, size = 2.5) +
theme_classic() +
theme_classic() +
xlab("Mean counts (log2)") +
ylab("Fold change log2(MANSC1/Mock)") +
ggtitle("MANSC1")
ggplotly(mansc_mock.maplot)Plot together
p = (snord_mock.maplot + ggtitle("") | mansc_mock.maplot + ggtitle("")) & ylim(-3, 0.75) & xlim(7, 11)
p
ggsave("figures/CD4-benchmark.pdf", p, width = 6, height = 4)MANSC1 vs. SNORD73
mansc_snord.res = compute_stats(ddsMat.fil, trt = "MANSC", ref = "SNORD", fitType = "parametric", lfcThreshold = 0.25, alpha = 0.05) %>%
left_join(annotations) %>%
mutate(log2FoldChange = round(log2FoldChange, 2)) %>%
mutate(log2FoldChangeShrunken = round(log2FoldChangeShrunken, 2)) %>%
mutate(gene_name = if_else(gene_name == "RPS3Amut", "SNORDmut", gene_name)) %>%
mutate(gene_name = if_else(gene_name == "RPS3Awt", "SNORDwt", gene_name))
# Get significant hits
mansc_snord.res %>%
filter(padj < 0.05)
mansc1_snord.res
# Browse
mansc_snord.res %>%
select(gene_id, log2FoldChange, padj, gene_name, gene_id) %>%
DT::datatable(extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
))
# Plot MA-plot
mansc_snord.maplot = mansc_snord.res %>%
ggplot(aes(x = log2(baseMean), y = log2FoldChange, label = gene_name)) +
geom_point(colour = alpha("grey", 0.7)) +
geom_point(colour = alpha("grey", 0.7)) +
geom_point(data = filter(mansc_snord.res, gene_name %in% c("MANSC1mut")), size = 3, colour = "#aec3a0") +
geom_point(data = filter(mansc_snord.res, gene_name %in% c("MANSC1wt")), size = 3, colour = "#a6a6a6") +
geom_text_repel(data = filter(mansc_snord.res, gene_name %in% c("MANSC1mut", "MANSC1wt")), box.padding = 0.5, max.overlaps = Inf, size = 2.5) +
theme_classic() +
theme_classic() +
xlab("Mean counts (log2)") +
ylab("Fold change log2(MANSC1/Mock)")
ggplotly(mansc_snord.maplot)Session info
sessionInfo()## R version 4.1.3 (2022-03-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices datasets utils methods
## [8] base
##
## other attached packages:
## [1] forcats_0.5.1 stringr_1.4.0
## [3] dplyr_1.0.9 purrr_0.3.4
## [5] readr_2.1.2 tidyr_1.2.0
## [7] tibble_3.1.7 tidyverse_1.3.1
## [9] ggpubr_0.4.0 plotly_4.10.0
## [11] ggplotify_0.1.0 readxl_1.4.0
## [13] DESeq2_1.34.0 SummarizedExperiment_1.24.0
## [15] Biobase_2.54.0 MatrixGenerics_1.6.0
## [17] matrixStats_0.62.0 GenomicRanges_1.46.1
## [19] GenomeInfoDb_1.30.1 IRanges_2.28.0
## [21] S4Vectors_0.32.4 BiocGenerics_0.40.0
## [23] ggrepel_0.9.1 ggplot2_3.3.6
## [25] ggsci_2.9 patchwork_1.1.1
## [27] knitr_1.39
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-3 ggsignif_0.6.3 ellipsis_0.3.2
## [4] XVector_0.34.0 fs_1.5.2 rstudioapi_0.13
## [7] farver_2.1.0 DT_0.25 bit64_4.0.5
## [10] AnnotationDbi_1.56.2 fansi_1.0.3 lubridate_1.8.0
## [13] xml2_1.3.3 splines_4.1.3 cachem_1.0.6
## [16] geneplotter_1.72.0 jsonlite_1.8.0 broom_0.8.0
## [19] annotate_1.72.0 dbplyr_2.1.1 png_0.1-7
## [22] BiocManager_1.30.18 compiler_4.1.3 httr_1.4.2
## [25] backports_1.4.1 assertthat_0.2.1 Matrix_1.4-0
## [28] fastmap_1.1.0 lazyeval_0.2.2 cli_3.3.0
## [31] htmltools_0.5.2 tools_4.1.3 gtable_0.3.0
## [34] glue_1.6.2 GenomeInfoDbData_1.2.7 Rcpp_1.0.8.3
## [37] carData_3.0-5 cellranger_1.1.0 jquerylib_0.1.4
## [40] vctrs_0.4.1 Biostrings_2.62.0 nlme_3.1-155
## [43] crosstalk_1.2.0 xfun_0.30 rvest_1.0.2
## [46] lifecycle_1.0.1 renv_0.15.4 rstatix_0.7.0
## [49] XML_3.99-0.10 zlibbioc_1.40.0 scales_1.2.0
## [52] hms_1.1.1 parallel_4.1.3 RColorBrewer_1.1-3
## [55] yaml_2.3.5 memoise_2.0.1 yulab.utils_0.0.5
## [58] sass_0.4.1 stringi_1.7.6 RSQLite_2.2.15
## [61] highr_0.9 genefilter_1.76.0 BiocParallel_1.28.3
## [64] rlang_1.0.2 pkgconfig_2.0.3 bitops_1.0-7
## [67] evaluate_0.15 lattice_0.20-45 labeling_0.4.2
## [70] htmlwidgets_1.5.4 bit_4.0.4 tidyselect_1.1.2
## [73] magrittr_2.0.3 R6_2.5.1 generics_0.1.2
## [76] DelayedArray_0.20.0 DBI_1.1.2 mgcv_1.8-39
## [79] pillar_1.7.0 haven_2.5.0 withr_2.5.0
## [82] survival_3.2-13 KEGGREST_1.34.0 abind_1.4-5
## [85] RCurl_1.98-1.6 modelr_0.1.8 crayon_1.5.1
## [88] car_3.0-13 utf8_1.2.2 tzdb_0.3.0
## [91] rmarkdown_2.14 locfit_1.5-9.6 grid_4.1.3
## [94] data.table_1.14.2 blob_1.2.3 reprex_2.0.1
## [97] digest_0.6.29 xtable_1.8-4 gridGraphics_0.5-1
## [100] munsell_0.5.0 viridisLite_0.4.0 bslib_0.3.1